library(hydroGOF)
library(clusterSim)

data=read.csv("C:/data/all_daily_merge.csv")
data=subset(data,data$station!=3)

data$aod=data.Normalization(data$aod,type="n5",normalization="column")
data$temp=data.Normalization(data$temp,type="n5",normalization="column")
data$avg_temp=data.Normalization(data$avg_temp,type="n5",normalization="column")
data$avg_rh=data.Normalization(data$avg_rh,type="n5",normalization="column")
data$avg_preci_24=data.Normalization(data$avg_preci_24,type="n5",normalization="column")


model=lm(pm25~aod+temp+avg_temp+avg_rh+avg_preci_24,data)
nSample=nrow(data)
nPredict= length(model$coefficients)-1

cutoff=4/(nSample-nPredict-1)
cook_thresold=qf(0.2,nPredict+1,nSample-nPredict-1)
cookValue=cooks.distance(model)
data["cookValue"] = cookValue

data1=data
data2=subset(data,data$cookValue<=cutoff)
data3=subset(data,data$cookValue<=cook_thresold)

for (i in 2010:2014) {
	testyear = i
	trainData1=subset(data1,data1$year!=testyear)
	testData1=subset(data1,data1$year==testyear)
		
	model1=lm(pm25~aod+temp+avg_temp+avg_rh+avg_preci_24,trainData1)
	pm25model=predict(model1,testData1)
	
	
	model_re=sum(abs(testData1$pm25-pm25model)/testData1$pm25)/nrow(testData1)*100
	model_rmse=rmse(testData1$pm25,pm25model)
	model_r=cor(testData1$pm25,pm25model)
	
	print(paste("Year model 1: ",testyear))
	print(paste("Number of training: ",nrow(trainData1)))
	print(paste("Number of testing: ",nrow(testData1)))
	print(paste("R2:",round(model_r*model_r,3)))
	print(paste("RMSE:",round(model_rmse,3)))
	print(paste("RE:",round(model_re,3)))

	
	print("------------------")
	
	
	trainData2=subset(data2,data2$year!=testyear)
	testData2=subset(data2,data2$year==testyear)
	
	model2=lm(pm25~aod+temp+avg_temp+avg_rh+avg_preci_24,trainData2)
	pm25model=predict(model2,testData2)
	
	
	model_re=sum(abs(testData2$pm25-pm25model)/testData2$pm25)/nrow(testData2)*100
	model_rmse=rmse(testData2$pm25,pm25model)
	model_r=cor(testData2$pm25,pm25model)
	
	print(paste("Year model 2: ",testyear))
	print(paste("Number of training: ",nrow(trainData2)))
	print(paste("Number of testing: ",nrow(testData2)))
	print(paste("R2:",round(model_r*model_r,3)))
	print(paste("RMSE:",round(model_rmse,3)))
	print(paste("RE:",round(model_re,3)))
	
	print("------------------")
	
	trainData3=subset(data3,data3$year!=testyear)
	testData3=subset(data3,data3$year==testyear)
	
	model3=lm(pm25~aod+temp+avg_temp+avg_rh+avg_preci_24,trainData3)
	pm25model=predict(model3,testData3)
	
	
	model_re=sum(abs(testData3$pm25-pm25model)/testData3$pm25)/nrow(testData3)*100
	model_rmse=rmse(testData3$pm25,pm25model)
	model_r=cor(testData3$pm25,pm25model)
	
	print(paste("Year model 3: ",testyear))
	print(paste("Number of training: ",nrow(trainData3)))
	print(paste("Number of testing: ",nrow(testData3)))
	print(paste("R2:",round(model_r*model_r,3)))
	print(paste("RMSE:",round(model_rmse,3)))
	print(paste("RE:",round(model_re,3)))
	print("------------------")
}

